module Distribution.PowerLaw
        ( powerlaw
        , estimatePower
        , PowerLaw
        )
        where

import qualified Numeric.Probability.Distribution as D

-- we made a significant assumption that x_min = 1 and we shift from [0..] to [1..]

shapePower k x = (x+1)**(-k)

powerlaw k = D.shape (shapePower k)

data PowerLaw = PowerLaw Double
        deriving (Show, Eq)

estimatePower :: [(Double, Double)]->Double
estimatePower xs        = (1.0+n/nsum)
        where  f (x0, y0) (x, y)        = (x0+x, y0+x*y)
               (n, nsum)                = foldl f (0.0, 0.0) xs

